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ABSTRACT 

In  Ballistic  Research  Laboratories  Memorandum  Report  No.  530, 

H Lotto  Method  of  Computing  Kill  Probability  of  Large  Warheads’1, 

F.  G.  King  describes  a random  sampling  method  for  determining  kill 
probabilities  of  a large  warhead  against  an  airplane.  To  obtain  the 
results  the  method  requires  a physical  model,  hand  drawing  of  random 
numbers,  and  the  use  of  kill  probability  curves  for  the  vulnerable 
components  of  the  airplane. 

In  this  report  a mathematical  model  for  the  purpose  of  solving 
the  problem  on  a high-speed  digital  computing  machine  is  presented. 
This  model  is  based  on  J.  von  Neumann's  suggestion  that  the  airplane 
be  replaced  by  several  ellipsoids  resembling  the  fuselage,  wings,  and 
engines.  The  necessary  formulas  for  computation  are  derived  from  the 
basic  geometric  model. 

The  kill  probabilities  are  determined  by  three-dimensional  inte- 
grals which  are  evaluated  either  by  random  sampling  methods  or  by 
strai^it forward  numerical  quadratures.  These  methods  are  compared 
from  the  viewpoints  of  accuracy,  speed,  and  machine  storage  require- 
ments. Limited  comparison  of  the  results  and  some  remarks  about 
applicability  of  more  general  problems  are  also  made. 

The  method  of  generation  of  the  pseudo-random  numbers  used  in 
the  random  sampling  procedures  is  also  described  in  the  appendix. 
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I.  INTRODUCTION 


1.1  Description  of  the  Problem 


An  airplane  may  be  brou^it  down  by  a bursting  warhead  by  damage  to 
its  structure  by  blast  or  by  destruction  by  fragments  of  the  warhead 
casing  of  a sufficient  number  of  vital  components , which  may  include 
pilots,  without  which  the  airplane  could  not  maintain  controlled  flight. 

When  the  criteria  for  a kill  for  a particular  airplane  are  estab- 
lished, using  basic  laws  for  combining  probabilities,  theoretically  one 
could  get  the  probability  of  obtaining  a kill  on  the  airplane  with  a 
particular  warhead  from  a knowledge  of  the  following: 


a)  the  joint  probability  distribution  of  the  burst  point  of 
the  warhead,  and  the  velocity  (speed  and  direction)  of  the  missile  at 
burst  point, 

b)  the  velocity  of  the  airplane  at  the  time  of  burst, 

c)  the  probability  of  destruction  of  the  structure  by  blast 
as  a function  of  directed  distance  of  the  missile  from  the  airplane, 
relative  velocities  of  the  missile  and  airplane,  wind  and  air  density 
at  time  of  burst, 

d)  the  static  distributions  of  fragment  size,  shape,  weight, 
and  velocity  at  time  of  burst, 


e)  the  aerodynamic  drag  on  the  fragments  as  a function  of 
the  size  and  shape  of  the  fragments  and  of  the  air  structure  at  the 
time  of  burst, 


f )  the  probability  distribution  of  the  air  structure, 
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g)  the  location  of  vital  components  and  of  shielding  com- 
ponents on  the  airplane,  and 


h)  the  probability  of  destruction  of  each  vital  component 
as  a function  of  direction,  size,  shape,  and  striking  energies  of 
fragments. 


Obviously,  for  practical  considerations,  such  as  the  impossibility 
of  getting  some  of  the  above  desired  probability  distributions,  the 
crudeness  of  those  obtained,  the  arbitrariness  of  the  definitions  of 
kills,  the  difficulty  of  computation  and  the  length  of  time  necessary 
to  carry  out  the  computations,  one  cannot  compute  the  probability  of 
such  a kill  exactly.  Therefore,  simplifying  physical  assumptions  were 
made  to  obtain  approximate  results.  This  is  justified  since  the 
available  vulnerability  data  for  this  problem  are  known  to  at  best  from 
ten  to  twenty  percent  and  the  errors  of  approximation  in  the  mathematic- 
al work  will  not  be  as  great. 
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Therefore,  in  the  particular  problem  discussed  in  this  report  it 
was  assumed  that: 

a)  The  missile  and  the  airplane  move  in  parallel,  horizontal 
planes  and  the  projection  of  the  direction  of  the  missile  on  the  plane 
of  motion  of  the  airplane  forms  a line  Ug°  off  the  nose  of  the  airplane. 

b)  The  distribution  of  burst  points  is  a trivariate  normal 
distribution  in  Cartesian  coordinates  (x,y,z)  with  a mean  (m-pm^m^) 

somewhere  in  the  plane  of  the  airplane  and  k5°  off  its  nose  in  the 

2 2 2 

direction  of  the  approach  of  the  missile  and  variances  <r^,  Ogs  ***■  a3 

and  zero  correlations.  (In  this  problem  the  origin  of  the  coordinate 
system  is  at  the  center  of  gravity  of  the  airplane ; the  x-axis  points 
to  the  starboard j the  y-axis  points  aft|  and  the  z-axis  points  aloft.) 

c)  The  velocity  of  the  missile  is  fixed  and  the  velocity  of 

the  airplane  is  negligible  relative  to  the  missile  and  fragment  veloc- 
itiesv 


d)  There  is  a region  called  the  "blast*1  region  about  the 
airplane  within  which  if  a burst  occurs  there  will  result,  with  prob- 
ability equal  to  one,  a kill.  Outside  this  region  the  probability  of  a 
kill  is  assumed  to  be  a function  solely  of  the  probabilities  of  des- 
troying individual  vital  components. 

e)  The  dynamic  distribution  of  lethal  fragments  is  confined 

to  a conical  region  about  the  nose  with  axis  coincident  with  the  missile's 
longitudinal  axis  and  to  a region  about  the  side  which  is  bounded  by 
cones  whose  axes  coincide  with  the  missile 8 s^longitudinal  axis.  Further- 
more, it  is  assumed  that  fragments  are  identical  in  mass,  size,  shape, 
and  in  resultant  speed  at  burst  and  travel  in  straight  lines  after  burst. 
It  is  also  assumed  that  the  dynamic  distribution  of  number  of  fragments 
per  solid  angle  subtended  from  the  burst  point  is  uniform  in  the  regions 
of  fragmentation  mentioned. 

We  remark  that  the  assumptions  a)  and  b)  are  not  restrictions  im- 
posed by  the  Lotto  method  or  by  any  of  the  methods  described  in  this 
report \ they  are  assumptions  of  the  particular  problem  solved.  The 
missile  could  move  in  any  other  straight  line  path  and  the  probability 
distribution  of  the  burst  points  could  be  different  without  any 
serious  modification  of  the  mathematical  formulation  for  the  ORDVAC  or 
of  the  coding. 

The  assumption  c)  is  not  as  crude  as  it  appears  in  its  statement, 
for  the  relative  motion  of  the  airplane  is  actually  considered  in  the 
problem  by  a modification  of  the  input  vulnerability  data.  • 
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1.2  The  Lotto  Method  of  Solution. 

A random  sampling  procedure  for  estimating  the  probability  of 
a kill  of  an  airplane  by  a warhead  under  the  assumptions  a)  througi  e) 
given  at  the  end  of  the  previous  section  has  been  successfully  used.  It 
is  described  in  reports  by  F.  G.  King  C2]  and  by  Stanley  Sacks  and  F.  G, 
King  [U]  from  which  some  of  the  introductory  material  in  this  report  is 
taken. 

Briefly,  the  procedure  is  as  follows.  A table,  called  a firing 
record  table  is  madej  the  headings  of  its  columns  include  "blast”,  the 
names  of  all  the  vital  components,  and  "kill".  A set  of  three  indepen- 
dent random  drawings  is  made  from  three  different  normal  populations 
whose  means  are  m,  , nu,  and  ra~,  respectively  and  whose  variances  are 
2 2 2 i-  <■  i 

aX3  02,  °3*  respectively.  This  triplet  determines  the  burst  point.  The 

axes  of  the  distributions  do  not  necessarily  coincide  with  those  of  the 
airplane  but,  in  general,  could  be  brought  into  coincidence  by  a ro- 
tation. Using  a scale  model  of  the  airplane,  it  is  first  determined 
whether  the  burst  is  sufficiently  close  to  the  airplane  to  destroy  its 
structure  by  blast*  If  so,  a one  is  tabulated  under  "blast" ; if  not,  a 
zero  is  scored.  Then  for  each  vital  component  successively  the  proba- 
bility of  destroying  it  is  determined.  This  probability  is  zero  if  the 
component  is  outside  the  lethal  fragmentation  region  of  the  missile  or 
if  the  component  is  shielded  from  the  lethal  fragment  spray  by  some  part 
of  the  airplane,  e.g.,  fuselage,  wings,  or  engines  ; otherwise,  this 
probability  is  assumed  to  be  a function  solely  of  the  distance  of  the 
component  from  the  burst  and  is  given  by  graphs  drawn  from  the  vulner-  .. 
ability  data  for  each  component.  Then  a random  number  from  a uniform 
population  in  (0,l)  is  drawn;  if  this  number  exceeds  the  above  proba- 
bility a zero  is  scored  in  the  firing  record  table  under  the  heading  of 
this  component;  otherwise,  a one  is  scared.  After  this  has  been  done  for 
all  components  a one  or  a zero  is  recorded  under  the  column  headed  "kill". 
A one  indicates  that  ones  appear  either  under  blast  or  under  a sufficient 
number  of  vital  components  of  the  same  type  which  if  destroyed  would  re- 
sult in  a kill;  a zero  indicates  the  negative  of  this.  „„ 

This  procedure  is  then  repeated  one  hundred  times  using  a new  ran- 
dom burst  point  each  time.  One  hundredth  of  the  total  number  of  ones 
in  the  kill  column  is  then  an  estimate,  in  the  sense  of  the  strong  law 
of  large  numbers  of  the  theory  of  probability  of  a kill  of  the  airplane 
by  "the  warhead.  This  estimate,  being  a random  variable  with  a binomial 
distribution,  has  a variance  given  by  pq/100  where  p + q ® 1 and  p is 
the  probability  of  a kill*  This  gives  a ratio  of  standard  deviation 

to  mean  of  0.1  -yj 4/p,  which  is  between  $%  and  20 % when  kill  probabilities 
is  between  0,8  and  0.2.  Thus  the  percentage  random  error  in  this  method 
of  computation  is  comparable  to,  actually  slightly  better  than,  the 
accuracy  of  the  data. 

Random  estimates  of  the  percentages  of  kills  due  to  blast  or  of 
kills  due  to  the  destruction  of  a sufficient  number  of  vital  components 
of  a given  type  as  well  as  estimates  of  the  probabilities  of  destroying 
ary  particular  component  or  combinations  of  components  are  also 
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obtainable  from  these  firing  record  tables.  Furthermore,  if  the  criteria 
on  numbers  of  vital  components  of  any  type  to  be  destroyed  to  get  a kill 
are  changed,  the  probability  of  a kill  is  obtainable  from  the  same  firing 
table.  In  addition,  using:  several  firing  record  tables,  one  is  able  to 
estimate  kill  and  other  probabilities  due  to  several  missiles  without 
constructing  new  tables.  In  fact,  the  estimates  of  these  probabilities 
could  be  obtained  from  one  firing  table  but  with  higher  variances j for  k 
missiles  the  variances  would  be  slightly  less  than  k times  the  variances 
for  one  missile.  This  is  so  because  instead  of  having  100  trials  in  the 
average,  there  would  only  be  100/k  trials  in  the  average.  The  presence 
of  multiply  vulnerable  components  actually  raises  the  kill  probability 
somewhat  and  this  lowers  the  expected  relative  error. 


H.  GRDVAC  SOLUTION  OF  THE  PROBLEM 


II. 3 The  Mathematical  Model. 


In  view  , of  the  simplifying  assumptions  a)  through  e)  of  1.2,  the 
mathematical  quantity  estimated  by  the  Lotto  method  is  the  integral 


(3.1) 


I » 


OO  oo 

^ di2(y)  f(x,y,z)  d#j(s) 


Pfrivtw) 


where  f(x,y,z)  is  the  probability  that  a kill  has  been  produced  (or  a 
blast  kill,  or  any  particular  combinat^m  of  vital  components  has  been 
destroyed)  by  a burst  at  (x,y,z)  and  dj^x)  dff  2(y)  <**3(2)  is  the 

probability  that  a burst  occurs  at  (x,y,z).  Ey  the  assumption  b)  of 


1.2,  we  have 
d^1(x)  d$2(y) 


(2n) 


1 

w. 


- 

e 


X“®1 


) 


& 


dx  dy  dz. 


1°2°3 


We  shall  estimate  the  same  integral  by  three  different  methods  which 
will  be  described  in  II.6.  Each  of  the  three  methods  has  in  ccomnon  the 
same  method  of  estimation  of  the  integrand  f(x,y,z)  which,  in  addition 
to  depending  on  the  burst  point,  (x,y,z),  depends  on  the  geometrical 
representation  of  the  airplane  and  the  vulnerability  data. 


Criteria  for  a good  geometrical  model  are  that 


a)  it  should  represent  the  airplane.  Its  more  bulky  components, 
and  the  blast  region  fairly  well, 

b)  It  should  contain  only  simple  formulas  that  do  not  require 
too  many  registers  in  the  machine’s  internal  storage  and  that  do  not  take 
much  computing  time  to  be  evaluated, 

c)  it  should  lead  to  simple  formulas  for  the  regions  shielded 
by  bulky  parts  of  the  airplane,  and 
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d)  it  should  lead  to  simple  determinations  of  whether  or  not 
a burst  is ' in  the  blast  region  or  a vital  component  is  in  a shielded 
region. 


Three  models  have  been  considered  as  approximations  to  the  form  of 
the  airplane  and  blast  regions  A set  of  ellipsoids  to  represent  the 
fuselage 5 wingss  engines,  and  blast  regionj  a set  of  spheres  to  represent 
the  above $ and  a set  of  cylinders  to  represent  the  above  except  the  wings 
which  could  be  represented  by  a rectangular  parallelipiped . They  all 
satisfy  the  above  four  criteria  fairly  well.  The  best,  however,  seems 
to  be  the  first,  idle  suggestion  of  J.  von  Neumann.  It  seems  to  satisfy 
a)  better  than  either  of  the  other  two.  All  three  models  satisfy  b) 
equally  well  except  that  many  more  spheres  than  ellipsoids  are  necessary 
to  represent  approximately  the  volumes  in'  question  and,  therefore,  require 
more  machine  storage  and  computing  time.  All  three  involve  only  quadric 
surfaces  which  require  only  a few  products  and  sums  to  describe  the  ele- 
mental surfaces  and  for  which  it  involves  only  simple  discriminations  on 
the  sign  of  a quadratic  or  a linear  form  to  determine  whether  or  not  a 
point  is  enclosed  by  the  surface.  To  determine  the  regions  shielded  by 
bulky  parts  of  the  airplane  is  more  difficult  using  the  third  model  than 
for  the  first  two.  Thus,  in  view  of  these  considerations,  the  ellipsoid 
model  was  chosen. 

Another  modification  is  the  replacement  of  volume-occupying 
vulnerable  components  by  point  components.  The  spray  regions  are 
^boundedTiy^cones  as  in  the  Lotto  method,  A vital  component  in  this 
region  is  considered  to  be  shielded  by  a part  of  the  airplane  represented 
by  an  ellipsoid  if  it  lies  outside  the  ellipsoid  and  lies  within  the 
region  inside  the  part  of  the  cone  determined  by  the  burst  point  and  the 
ellipsoid  which  is  on  the  other  side  of  the  ellipsoid  from  the  burst 
point.  The  mathematical  criteria  for  this  will  be  derived  in  the  next 
section. 

II. U Derivation  of  Mathematical  Criteria  for  Determining  Whether  or  Not 

a Vital  Component  is  in  the  Lethal  Fragnent  Spray. 

The  probability  that  a vital  domponent  at  a point  C is  destroyed 
by  fragnents  from  a burst  at  a point  P is  zero  if  the  component  is  not 
in  the  lethal  fragnent  spray  and  is  a function  of  the  distance  |C  - P| 
if  the  component  is  in  the  lethal  fragnent  spray.  A vital  component  at 
C (See  Fig.  l)  is  assumed  to  be  vulnerable  to  fragments  from  a burst  at 
P of  a missile  whose  unit  direction  vector  is  V at  the  time  of  burst  if 
it  is  either  inside  the  nose  spray  cone  or  outside  both  and  the 

rear  cone  Kg  and,  in  addition  to  either  of  these  conditions,  is  not 

shielded  or  masked  from  the  burst  by  some  part  of  the  airplane  repre- 
sented by  ellipsoids,  one  of  which  is  represented  in  Fig.  1 by  E.  If  C 
is  inside  E it  is  still  considered  vulnerable  if  it  is  in  the  spray 
regions  and  not  shielded  by  another  ellipsoid.  (Perhaps  a better 
assumption  would  be  that  if  C is  inside  E it  is  vulnerable  only  if -it 
is  on  the  same  side  of  the  polar  plane  of  P with  respect  to  E as  P. 

" However,  the  way  in  which  the  data  for  the  kill  probabilities  as  a 
function  of  distance  are  gathered  and  averaged  militates  against  this.) 
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FIGURE  1 


Fig.  1 represents  a dynamic  diagram  and,  therefore,  the  angle  p is 
greater  than  the  angle  a because  the  momentum  of  the  missile  is  added  to 
the  momentum  of  the  fragments  acquired  in  the  burst. 

To  obtain  the  mathematical  criteria  for  determining  whether  a 
vital  component  is  or  is  not  in  the  lethal  fragnent  spray  we  need  the 
stations  for  the  cones  Kp  Kg,  and  (K^  which  we  shall  derive  directly 

from  vector  considerations.  Z 

A point  X ® (xp  Xg,  x^)  is  on  either  nappe  of  the  cone  Kp  one  of 
whose  nappes  is  represented  in  Fig.  1 if 

^ v 

-^(x)  - [(X  - P)  . v]2  - |x  - P I2  cos2  0 = 0. 

It  is  outside  the  nappes  of  the  cone  if  the  quadratic  form  K.  (X)  is 
positive  and  it  is  inside  if  K^(X)  is  negative.  Similarly,  a point  X 
is  on  either  nappe  of  the  cone  Kg,  is  outside  the  nappes  of  Kg,  or  is 
inside  the  nappes  of  Kg  according  as  the  qiadratic  form 

-Kg(X)  - [(X  - P)  • V]2  - |X  - P|2  cos2  p 

is  zero,  is  positive,  or  is  negative,  respectively. 

A point  Y ji  P is  on  the  cone  K^)  the  cone  determined  by  the  pencil 

of  lines  through  the.._bur.st„point  P and  tangent  to  the  shielding  ellip- 
Spid  'E,  if  and  only  if  the  line  through  P and  Y is  one  of  these  tangent 
lines.  To  express  this  mathematically,  let  us  assume  for  the  while  that 
R the  center  of  the  ellipsoid  E is  also  the  origin  of  our  coordinate 
system,  and  let  X *=  (x^,  Xg,  x^)  be  a point  on  the  ellipsoid  E whose 

eqaation  in  this  system  is 

(U.2)  XAXT  ® 1 

T 

where  X is  the  transpose  of  the  row  vector  X and 
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In  which  a,  b,  and  c are  the  principal  semi-axes  of  the  ellipsoid  E. 

If  I is  a point  on  E and  t is  a parameter , then 

(U-3)  Y - P (l~t)  + t X 

is  the  equation  of  a line  through  P and  X®  Now  Y is  on  the  ellipsoid 
E if 

[p(l-t)  + t x]a  [p(l-t)  + t x]T  - 1 

or,  since  A is  symmetric, 

Oi.U)  t2(PAPT  - 2 PAXT  + XAXT)  ♦ 2t(PAXT  - PAPT)  ♦ (PAP1  - l)  » 0. 

Now  I is  a tangent  point  on  E if  and  only  if  there  is  a unique  t 
satisfying  (U.U),  i.e*,  if 

(paxt  - papt)2  - (PAPT  - 2 PAXT  + XAXT)  (PAPT  -i  1)  = 0, 
which  when  expanded  and  rearranged  yields 

(U.5)  (paxt  - i)2  - (papt  - 1)  (XAXT  - 1)  - 0® 

But  since  X is  on  E (U«£)  reduces  to 
(U«6)  PAXT  - 1. 


This  means  that  X must  lie  on  the  polar  plane  of  P with  respect  to  E® 

To  obtain  the  eqiation  of  the  cone  (Obviously,  P must  be  outside  E 

for  to  be  a real  cone®),  we  let  the  parameter  t in  (lu3)  range  over 

the  real  numbers  and  the  variable  point  X range  over  the  ellipse  deter- 
mined by  the  ellipsoid  E and  the  polar  plane  of  P with  respect  to  E. 

To  determine  the  equation  of  this  cone  in  non-par ametric  form  we  elimin- 
ate X and  t between  (U»2),  (U*3),  and  (U»6)  for  t / 0,  since  t « 0 yields 
Y - P. 

Letting  t / 0 and  solving  Qi«3)  Yor  X,  we  obtain 
(U.7)  X - P + i (Y  - P) 

which  we  substitute  in  (I*°2 ) to  get 

(u«8)  PAPT  + | PA  (Y  - P)T  + (Y  - P)  A (Y  - P)T  - 1® 
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Substituting  (It  ,7.)  in  (it *6),  we  obtain 

(It .9)  PAPT  + - PA  (Y  - P)T  - 1 

Subtracting  (it* 9)  from  (I4.8)  and  multiplying  the  result  by  t,  we  have 

(lt.10)  PA  (Y  - P)T  + ^ (Y  - P)  A (Y  - P)T  - 0. 

Eliminating  l/t  between  (it»9)  and  (UolO),  we  arrive  at 

PAPT  (Y  - P)  A (Y  - P)T  - [PA  (Y  - P)^2  - (Y  - P)  A (Y  - P)T  = 0. 
Expanding  this  and  using  the  symmetry  of  A,  we  have 

papt  yayt  - 2 papt  payt  + (faft)2 

- (payt)2  + 2 payt  papt  - (papt)2 

- YAYT  + 2 PAYT  - PAFT  = 0 


or 


(PAPT  - 1)  (YAYT  - 1)  - (PAYT  - 1) 


2 


0 


for  the  equation  of  the  cone  Ky 

Thus  a point  X is  on  either  nappe  of /K^jis  outside  the  nappes  of  ■ 
or  is  inside  the  nappes  of  accordingas  the  quadratic  form 


(lt.ll) 


k3(x) 


(PAP1  - 1)  (XAXJ 


1)  - (PAX1  - if 


i 

is  zero,  is  positive,  or  is  negative  respectively.  For  each  shielding 
ellipsoid  there  is  a different  matrix  Aj  therefore,  we  have  a set  of 
quadratic  forms  of  the  type  of  K^.  If  there  are  K shielding  ellipsoids, 

E^,  let  us  call  their  respective  matrices  A^,  k ■ 1,  2,  ...,  K,  and  the 

quadratic  forms  associated  with  the  cone  determined  by  P and  the  k-th 
ellipsoid  will  be  designated  by  K*  , . Thus 

(14.12)  K^k00  ■ (PAjPT  - 1)  (XA^  - 1)  - (PA^X1  - 1)  . 


It  should  be  recalled  that  the  origin  of  the  coordinate  system  in  the 
formula  is  the  center  of  Ek  and  that  the  qiantities  P and  X are 

usually  measured  relative  to  the  center  of  gravity  of  the  airplane. 
Therefore,  the  appropriate  translations  should  be  made  before  using  the 
formula  (£.12 ).  We  return  now  to  our  original  reference  coordinate 
system  with  its  origin  at  the  center  of  gravity  of  the  airplane. 
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The  fragment  spray  kill  probability  for  each  component  is  assumed 
to  be  different  in  the  nose-spray  region  (inside  K^)  from  that  in  the 

side-spray  region  (outside  and  Kj)  but  is  assumed  to  be  uniform  with 

respect  to  the  angle  subtended  from  V (See  Fig.  1.)  in  each  region..  If 
the  point  representing  a particular  component  falls  in  the  nose-spray 
region  and  is  not  shielded  by  ary  ellipsoid  E^,  let  ue  call  the  procedure 

of  computing  the  probability  of  destruction  of  the  component  by  frag- 
ments procedure  A?  if  the  component  falls  in  the  side  spray  region  and 
is  unshielded,  the  corresponding  procedure  will  be  called  procedure  B? 
and  if  the  component  falls  inside  the  rear  cone.  Kg  or  is  shielded,  the 

corresponding  procedure  (which  is  simply  recording  a zero)  will  be 
called  procedure  C. 

The  logical  steps  to  determine  the  choice  of  the  appropriate  one  of 
the  above  procedures  are  as  follows.  First  the  scalar  product  (C  - P)  ® V 
is  formed . If  this  is  positive  then  the  component  C lies  on  the  same 
side  of  the  equatorial  plane  of  the  missile  at  burst  as  the  direction 
vector  of  the  missile  points ? if(C-P)  • T is  negative  then  C is  on  the 
other  side  of  the  equatorial  plane „ 

If  (C  - P)  » 7 > 0,  then  the  component  C is  either  in  the  nose 
spray  region  or  side  spray  region,  and  it  is  next  determined  whether  it 
is  shielded  by  one  or  more  of  the  ellipsoids  E^,  the  logic  of  which 

determination  will  be  outlined  later.  If  the  component  is  shielded,  we 
pass  on  to  apply  procedure  C.  If  it  is  not  shielded,  then  K^(C)  is 

formed.  If  K^(C)<  0,  procedure  A is  applied  since  C is  then  in  the 

nose  spray  and  is  vulnerable.  If  IC^(C)  0,  then  G is  in  the  side 

spray  region  and  procedure  B is  to  be  followed. 

If  (C-P)  *7^0,  then  C Is  either  in  the  side  spray  region  or 
in  the  rear  cone  region  of  no  fragment  spray.  In  order  to  determine  in 
which  of  these  two  regions  C lies,  Kg(C)  is  formed?  if  Kg(C)<  0,  C is 

in  the  no  spray  region  and  procedure  C is  applied?  if  Kg(C)  2 0,  it  is 

next  determined  whether  C is  shielded  or  not.  If  C is  shielded,  pro- 
cedure C is  then  applicable,  and  if  C is  unshielded,  we  go  on  to  apply 
procedure  B. 

When  C is  in  a spray  region,  in  order  to  determine  whether  it  is 
shielded  or  not  the  subsequent  logical  procedure  is  followed.  Re- 
calling that  the  center  of  the  k-th  shielding  elliposid  is 

k = 1,  2,  K,  we  form  (P  - \)  Afc(C  - Rk)T  - 1.  If 

(P  - R^)  Ajf(C  - R^)^  - 1<  0,  then  P and  C are  on  opposite  sides  of' 

the  polar  plane  of  P with  respect  to  the  ellipsoid  E^  and  we  fora 

K3,k{c>  - [<p-v  V^-V1  -1]  - [<p-v  vc-vT  -1!- 
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If  lc(C)<  0,  C is  inside  the  cone  k and  we  form 
(C-R^  Afc  (C-R^1  - 1.  If  (C-R^)  Aj^  - 1 > 0,  then  C lies  out- 

side the  ellipsoid  but,  as  a result  of  previous  discriminations,  it 
is  found  to  be  inside  k and  on  the  opposite  side  from  P of  the  polar 

plane  of  P with  respect  to  E,  • Hence  C is  shielded  and  we  pass  on  to  the 
procedure  C. 

On  the  other  hand,  if 
Afc  (C-Rfc)T  - 1 < 0 m 

0,  and  (C-Rk)  Aj£  (C-R^)7  -1^0,  then  C is  not  shielded  by  Ek 

because  it  will  be  on  the  same  side  as  P of  the  polar  plane  of  P with 
respect  to  or  it  will  be  on  the  other  side  but  also  outside  the  cone 

K,  , or  it  will  be  inside  K.  . but  also  inside  E,  and  thus  not  shielded 

according  to  our  assumption.  If  C is  then  not  shielded  by  Ek  and  k < K, 

k is  replaced  by  k + 1 and  the  procedure  is  repeated.  If  k « K then 
procedure  A or  B is  applied  depending  on  whether  C is  in  the  nose  spray 
region  or  the  side  spray  region. - 


(P-Rfc)  Ak  (C-Ty1  - lfeo,  or  if 
i K-  k(C)  26  0,  or  if  (P-Rj^)  Afc  (C-R^1  - 0, 


It  should  be  remarked  that  the  restriction  that  the  matrices  Ak  in 

(U.12 ) be  diagonal  is  particular  to  the  actual  problems  run  on  CRDVAC 
but  is  not  necessary  to  the  derivation  of  (U.12).  If  the  shielding 
ellipsoids  had  axes  which  were  not  parallel  to  the  coordinate  axes 

then  the  Ak  would  not  be  diagonal  but  (U.12)  would  still  be  valid.  For 
airplanes  with  swept  back  wings  such  modifications  are  necessary.  It 
is  possible  that  for  delta  wings  ellipsoidal  approximations  are  too 
crude.  In  which  case  a triangular  prism  of  very  small  depth  cculd  be 
used  with  a small  increase  in  memory  requirements  for  the  necessary 
discriminations  to  determine  the  shielded  regions  for  each  burst  point. 


The  logical  flow  chart  of  the  application  of  the  criteria  of  this 
section  for  one  component  is  gi^n^in-Eig.  2 . 

II. 5 Computation  of  the  (Kill  Pr ob ab ility^f^r  a Given  Burst  Point. 

The  integrand  f (x,y,zj  of" ( 3.1),  which  is  the  probability  that  a 
kill  has  been  produced  by  a burst  at  (x,y,z),  was  programmed  for  com- 
putation by  two  different  methods.  One  method  has  the  advantage  of 
being  very  general,  but,  on  the  hand,  with  our  specific  kill  criteria 
it  required  about  seven  times  the  computing  time  of  the  other.  The 
recent  addition  of  a new  order,  logical  "and11  to  the  CRDVAC1  s list  of 
instructions  will  by  its  use  make  the  two  computing  times . nearly  the 
same.  Both  methods  will  be  described  here. 


The  more  general  method  will  be  discussed  first.  If  there  are  M 
vulnerable  components  of  the  airplane,  then  when  a missile  explodes. 


IS 


M 

there  are  2 mutually  exclusive  possible  events,  besides  blast , corres- 
ponding to  the  destruction  or  non-destruction  of  each  of  the  components 
Ci,  Cg,  «9o,  Cjjo  Bach  of  the  possible  events  may  be  represented  by  an 

M . 

M-digit  binary  integer  £ 21  p*  where 


if  C±  is  destroyed 
if  C.  is  not  destroyed. 


Let  us  define  a partial  ordering  relation,  ( ^ ),  among  the  M 
binary  integers  such  that 


M . M 

X 21  p.  (*)  X 


i-1 


i°l 


21  y„ 


if  and  only  if  p^  < y^  for  each  i = 1,  2,  *99,  M.  Let  S be  the  subset 

of  the  2^  possible  M-digit  binary  integers  which  represent  kills . Let 

S be  that  subset  of  S such  that  a)  if  s is  in  5 then  there  exists  an 

element  s in  S such  that  s (S  ) s and  b)  if  s is  in  S then  there  does 

not  exist  another  element,  s,  of  S such  that  s(S)  s 9 If  r is  an  M- 
digit  binary  integer,  let  us  define  a function  € (r)  such  that 


<(r) 


r * 

l 0 1 


if  there  exists  an  s in  S such  that  s (-fi)r 


if  there  does  not  exist  an  s*  in  S such  that  s*(<)r. 


We  remark  that  fortf(r)  = 1 it  is  necessary  and  sufficient  that  r be  in 
S. 


The  probability  p^,  of  destruction  of  the  component  Ch  is,  in  the 

problem  handled,  a function  of  the  distance,  D,  frcrn  the  component  to 

the  given  burst  point  (x,y,z)  the  density  of  the  lethal  fragment  spray, 
and  the  vulnerable  area  of  the  component,— A— 9-~_It  is  given  by 

Pl  "I*®  / 

where  a .is  a constant  proportional  4x?-tfie^fragnent  density.  The  vulner- 
able area  is  obtained  from  e xper  imental~cfata  and'for  IffiDVAC  computations 

was  fitted  by  a curve  composed  of  segments  of  parabolas  and  straight  lines. 
This  fitting  was  done  by  W.  Barkley  Frrcz. 
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If  is  the  probability  of  destruction  of  the  canponent  by  a 

. i jf.-~  M 

burst  at  the  given  point  (x,y,z)  ana  if  r 1=  21  p,  (r),  then  the 

i-1 

probability,  P(r),  of  occurrence  of  the  event  represented  by  the  binary 
integer,  r,  is  given  by 

pm  - Pi + f1  - Mr)]  & - pi>} 


and  the  probability,  of  a kill  due  to  fragnents  from  a burst  at 
^x>y>z)  is  giran  by 

*-1 

PF  ■ S *(r)  P(r). 
r = 0 


This  method  is  very  general,  but,  for  the  specific  kill  criteria 
in  the  actual  problem  solved  on  QRDVAC,  a second  method  was  also  used  in 
which  it  was  possible  to  decrease  the  large  number  of  operations  re- 
quired by  avoiding  the  comparisons  involved  in  the  relation  (■<?)  and 
reducing  the  total  number  of  probabilities  needed  to  compute  Pp®  The 

efficiency  of  the  second  method  depends  on  the  fact  that  in  the  cases 
treated  the  set  of  vital  components  was  composed  of  a number  of  dis- 
joint subsets  containing  only  components  of  the  same  type  and  the  number 
of  any  one  type  was  smallj  furthermore,  for  each  type  there  exists  a 
fixed  number  of  vital  components  which  must  be  destroyed  in  order  to 

achieve  a kill.  If  there  are  J different  types  and  if  Q . is  the 

J 

probability  of  killing  at  least  the  required  number  of  vital  components 
of  the  type,  J^,  sufficient  to  produce  a kill,,  then 


T 

Pp  « 1 - TT  (1  - «,). 

J-l  . J 

Given  the  probabilities^  p^,  since  the  number  of  components  of  each 
type  is  small,  the  computation  Q.  is  very  simple.  For  example,  if  J- 
contained  only  the  components  and  C^,  and  if  the  destruction  of 
only  one  were  sufficient  for  a kill  then  % “ + Pg  “ P1P2*  ^ ^1 

contained  only  the  components  C^,  Cg,  and  C^,  and  if  the  destruction 

of  only  two  were  sufficient  for  a kill  then 
% - P-lP2  + PXP3  + P2P3  - 2p1p2p3* 

Finally,  if  is  the  probability  that  the  airplane^  structure 
is  destroyed  by  blast  from  a burst  at  the  given  point  (x,y,z),  then 


\ 
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the  integrand  of  (3«l)  is  given  by 


fl  if  (x,y,z)  is  in  the  blast  region 

Pp  if  (x,y,z)  is  not  in  the  blast 

- 'f  / region. 

11,6  Choices  of  Methods  of  Evaluation  of  the  Kill  Probability  Integral. 

Several  methods  of  quadrature  to  evaluate  the  kill  probability 
integrals  given  by  (3.l)  were  considered.  These  methods  fall  into  two 
general  classes,  random  sampling  or  Monte  Carlo  procedures  and  systematic 
numerical  integration. 

Criteria  for  choosing  any  method  over  others  are  usually  based  on 
the  accuracy  desired  in  the  value  of  the  integral  (3.1),  on  "the  accuracy 
and  speed  available  in  the  machine,  and  machine  internal  storage  require- 
ments for  each  method.  In  this  problem  since  sufficiently  accurate  and 
sufficiently  fast  methods  with  about  the  same  machine  internal  storage 
requirements  are  available  in  both  the  above  two  general  classes, 
storage  requirements  are  of  minor  importance  in  the  making  of  choices 
among  the  methods  considered.  The  programs  actually  used  occupied 
nearly  all  of  the  GRDVAC  storage,  but,  if  it  were  necessary,  the  same 
procedures  could  be  coded  to  take  a somewhat  smaller  storage  - about 
850  or  so  words. 

In  random  sampling  procedures  since  a correct  statement  of  the 
answer  is  that  the  answer  is,  say  x with  probability  p,  the  criterion 
assigning  greater  accuracy  to  a particular  random  sampling  method  than 
to  another  usually  that  the  first  method  have  a smaller  variance  than 
the  second.  In  choosing  one  sampling  method  over  another  we  shall  use 
this  criterion  in  addition  to  the  criterion  of  speed. 

Among  the  random  sampling  integrating  procedures  three  were  given 
some  consideration  for  possible  use  on  GRDVAC;  let  them  be  called 
RSIP-I,  RSIP-II,  and  RSIP-III. 

RSIP-I  is  the  Lotto  procedure  described  in  1,2  with  the  mathe- 
matical model  given  in  11,3. 

3h  RSIP-II  a seqience  of  N points  { (xn,  yn,  zQ)  n * 1,  2,  N, 
is  chosen  from  a trivariate  normal  distribution  with  mean  (nip  m^,  m^)  and 
variance  (crp  erg,  cr^).  Then  an  estimate,  in  the  sense  of  the  strong  law 
of  large  numbers,  of  the  integral  (3,1)  is  given  by 
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In  RSEP-III  a sequence  of  N points  £ (xr,  yn,  zn^  n = 1,  2,  N, 

is  chosen  from  a trivariate  uniform  distribution  in  the  region  defined  by 

< ka-L»  | y-™2  I - |a“m3| 


Then  an  estimate,  in  the  sense  of  the  strong  law  of  large  numbers,  of  the 
integral  (3.1)  is  given  by 


The  ranges  of  8op  fJc^,  and  6a  ^ far  the  random  variables  x^,  y..,  and  z^ 

were  chosen  because  the  probability  mass  within  iio  of  the  mean  in  the 
univariate  normal  distribution  is  0.9999365*  Since  accuracies  of 
several  percent  in  the  answers  are  more  than  are  necessary  in  this 
problem,  ranges  of  seven  a*s  or  even  six  a's  could  be  used  instead  of 
eight  a's*  In  the  normal  distribution  the  probability  mass  within  3.5c 
of  the  mean  is  0.999535  and  it  is  0.997300  within  3a  of  the  mean. 

The  random  number  generation  used  in  these  schemes  is  discussed  in 
the  Appendix. 

Among  the  systematic  numerical  integrating  procedures  two  were 
considered;  let  them  be  called  SNIP-I  and  SNIP-II.  Because  the 
accuracy  desired  in  the  problem  is  of  the  order  of  several  percent, 
simple  procedures  are  sufficient.  Therefore,  each  of  these  two  pro- 
cedures considered  were  three-dimensional  Riemann  sums. 

In  SNIP-I  the  summands  are  evaluated  at  Ir  points  at  the  geometric 
centers  of  elemental  cubes  of  equal  volume  in  the  physical  space.  For 
the  same  reasons  expressed  in  the  description  of  RSIP-III  the  region  of 
integration  is  restricted  to 


I*-™]. I - rol»  |y^2|  ^ **29  | z-m3  | * 
where  r may  be  chosen  as  3,  3*5>  or  U*  The  integral  (3.1)  is  estimated 


by 


*»m_  2 


4 0 p p p e 2 °1 

T i»l  .1=1  k=l  d 


. i * (2a, 


♦ (■ 


where 


/ ' 2i-lv 
m,  - ro.  U - -» — ) 


7^  = mg  - ra2  (l  - 
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and 


•k  ' m3  ‘ “3  (1  ' T-*- 

In  SNIP-II  the  summands  are  evaluated  at  L points  at  the  proba- 
bility mass  centers  of  elemental  cubes  of  equal  probability  mass  in  a 
three-dimensional  probability  space  whose  probability  mass  distribution 
is  a trivariate  normal  distribution.  The  integral  (3.1)  is  estimated  by 


where 


T t E f 

■ ■ ■ >1  ■ 


i"l 


k=l 


fCn^  + cr1ti,  mg  + Ojjtj,  “3  + <*3*^) 


ami  for  v £s  L/2 


o 


II. 7 Comparison  of  Methods  of  Evaluation  of  the  Kill  Probability 
Integral. 

In  this  section  the  methods  listed  in  II. 6 will  be  compared  on 
the  basis  of  theoretical  accuracy,  speed  and  memory  requirements. 

The  random  sampling  integrating  procedures  will  be  considered 

first. 

If  {x/^  is  a sequence  of  independent  identically  distributed 

random  variables  and  G(x)  is  a function  of  a very  general  class  (Baire 
functions  are  included  in  this  class.),  then 

Varfn  S Q(Xi>}  “KVar  {Q(Xi>}  "nE  f^i*]  “n®2  f °<Xi>  ^ . 

This  formula  will  be  used  to  determine  the  respective  theoretical 
variances  of  the  three  random  sampling  integrating  procedures.  The  re- 
spective variances  will  be  designated  as  Varj,  Varjp  and  Varjjj. 

In  RSIP-I,  G(X^)  is  a Bernoulli  random  variable  which  takes  on  the 

value  0 with  probability  I and  the  value  1 with  probability  1-1,  where 
I is  the  kill  probability  integral  (3.1).  Thus 

(7.13)  Tarz  = | (I-I2). 
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In  RSIP-II,  Xi  is  a three-dimensional  random  variable  with  a tri- 
variate normal  distribution  witii  mean  ecjial  to  (m^,  mg,  m^)  and  variances 
equal  to  (o^,  Cg,  o^)  and  G^)  - f(x^,  j^,  zi),  the  integrand  in  (3*1). 
Therefore, 


’ oo  oO 


(7*14)  Varn  “ ^ | \ d^Cx)  \ d $ g (y ) t ^(x^z)  d^(z)  * I2 

- C>0  «*06 


In  RSIP-III,  is  a three  dimensional  random  variable  with  a tri- 
variate uniform  distribution  in  .... 


jx-m-jJ  < rap  (y-nigl  ^rog,  <,  ro^ 


where  the  value  of  r could  conveniently  be  3,  3*5,  or  4. 


Then 


G(Xi)  = (“0  f(xi,yi,zi)  e 


3/2 


5 [<¥»>■  • ep>‘  • i 
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and 
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o 2 3/2 
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.15)  VarII;r  = j ^ dfr^x)  ^ d$g(x)  ^ (— 

1 “ro2  **3 


(^)2 ♦ (^)2J  , K o 

a2  a3  f^x^z)  d$^(y)  - I j 


The  only  difference  in  the  expressions  for  VaTj  and  VaTjj  is  in  the 

power  of  the  integrand,  f(x,y,z),  in  the  first  members  of  the  right  hand 
sides  of  (7*13)  and  (7.14).  In  the  former  f(x,y,z)  appears  to  the  first 
power,  while  in  the  latter  it  is  squared.  Hence,  since  f(x,y,z)  is  the 
probability  of  a kill  or  some  other  desired  event  from  a burst  at  the 
point  (x,y,z)  and,  therefore,  is  between  zero  and  one,  Var S VaTj. 


In  the  Lotto  method  it  has  been  customary  to  choose  M « 100.  For 
this  the  root  mean  square  relative  error  is  0.1  y(l-IJ/I,  implying  that 
if  a$  is  a bound  on  the  allowable  root  mean  square  relative  error  -then 

(7*16)  (1  + flVlOO)-1  * I. 


This  implies,  for  example,  that  if  the  root  mean  square  relative  errors 
are  not  to  exceed  10$  then  I must  not  be  less  than  1/2;  if  the  root 
mean  square  errors  are  not  to  exceed  20$  then  I must  exceed  1/5.  (Of 
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course,  choosing  M to  be  larger  will  reduce  the  error,  so  that  if  small 
values  of  I are  to  be  computed  with  accuracy  more  extensive  sampling  is 
needed.)  To  test  the  random  numbers  used  as  well  as  to  compare  RSIP-I 
and  RSIP-II  fifty  runs  of  the  problem  were  made  using  RSIP-H  with 
different  and  independently  randomly  chosen  values  of  rQ  at  the  request 

of  F.  G.  King.  The  values  of  the  variances  obtained  for  almost  all  of 
the  probabilities  of  kills  of  component  combinations  indicated  that 
approximately  l.$  to  3*2  times  as  many  burst  points  would  have  been  re- 
quired to  produce  the  same  accuracy  if  RSIP-I  had  been  used  instead  of 
RSIP-II.  For  the  probabilities  of  kills  due  to  destruction  of  a 
sufficient  number  of  one  particular  type  of  vital  component  from  $ to  8 
times  the  number  of  burst  pointe  used  in  RSIP-II  would  have  been  needed  to 
get  the  same  accuracy  if  one  used  RSIP-I.  However,  since  these  proba- 
bilities are  small  and  contribute  only  a little  to  the  total  probability 
of  a kill,  great  accuracy  in  these  small  probabilities  is  certainly  not 
necessary  unless  one  were  strongly  interested  in  the  particular  com- 
ponents1 probability  of  destruction. 


The  difference  VarT  - VarTT  is  given  from  formulas  (7.13)  and 

(7.U.)  as  11 


Var j - Var^j 


| T dS^x)  J d$2(y)  { f(l-f)  d*2<«). 

■ eo  _CK3  .40 


Therefore,  from  this  formula  one  sees  that  if  most  of  the  probability 
mass  is  near  where  f(l-f)  is  near  its  maximum,  i.e.,  where  f = 1/2 , 
then  the  difference  Var^  - VaXjj  is  of  the  order  of  Var^.  In  the  case 

of  singly  vulnerable  components  or  multiply-vulnerable  components 
physically  close  together  the  bulk  of  the  total  kill  probability  on  a 
component  combination  is  obtained  from  points  near  the  components,  or 
near  the  component  combination  as  the  case  may  be,  where  f is  near  .1. 

The  storage  requirements  and  the  computing  time  required  for 
RSIP-I  is  about  the  same  for  RSIP-II  for  the  same  number  of  bursts 
because  of  the  extra  randomization  in  RSIP-I  while  their  printing 
times  are  identical  provided  that  the  sequences  of  zeros  and  ones  which 
are  the  outcomes  of  each  burst  in  RSIP-I  are  not  printed.  On  the  other 
hand,  if  these  zeros  and  ones  are  to  be  recorded  as  in'  the  hand  Lotto 
method  the  printing  time  in  RSIP-I  is  greater  than  in  RSIP-II.  Further- 
more, if  the  same  root  mean  square  relative  error  is  desired  in  RSIP-I 
as  in  RSIP-II,  as  has  been  pointed  out,  for  moderate  values  of  (3«l) 
at  least  twice  the  computing  time  of  RSIP-II  is  required.  Therefore, 
on  the  basis  of  these  considerations  RSIP-II  is  recommended  by  the 
present  authors  over  RSIP-I,  provided  that  firing  record  tables  are 
not  required  to  be  printed  for  use  in  other  problems. 

The  storage  requirements  for  RSIP-II  are  slightly  higier  than  for 
RS IP-in.  The  canputing  time  required  for  RSIP-IU  is  approximately  2/3 
that  of  RSIP-II  for  the  same  number  of  bursts.  This  is  caused  by  the 
fact  that  RSIP-II  requires  random  numbers  chosen  from  normal  distri- 
butions while  RSIP-III  requires  random  numbers  chosen  from  a uniform 
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distribution  and  the  method  used  on  ORDVAC  to  produce  normally  distri- 
buted random  numbers  transformed  uniformly  distributed  random  numbers 
into  normally  distributed  random  numbers  and  used  a routine  to  compute 
logarithms  as  part  of  the  transformation;  this  routine  requires  a good 
fraction  of  the  computing  time. 


However,  Var^  is  somewhat  less  than  Varjjj.  The  only  difference 

between  the  first  members  of  the  right  hand  side  of  (7*lli)  and  of  (7*15) 
is  that  in  (7*15)  there  appears  the  extra  factor 


(7.17) 
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this  factor  exceeds  unity  and  is  less  than  unity  elsewhere.  Furthermore, 

f2(x,y,z)  is  generally  larger  in  the  region  is  determined  by 
(7.18)  than  elsewhere  where  it  is  probably  very  small  because  it  is  the 
square  of  the  kill  probability  for  a burst  at  (x,y,z)  and  generally  de- 
creases away  from  the  origin  which  is  not  far  from  the  mean  burst  point, 
certainly  very  much  closer  than  5»2k  min  (opOg,^).  Finally,  since 
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and  because  of  the  foregoing  remarks  about  f2 (x,y,z ),  a large  fraction 
of  the  value  of  the  first  members  of  the  right  hand  sides  of  (7«lU)  and 
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(7*15)  is  contributed  by  the  integration  over  the  region  determined  by 
(7*18)  in  the  case  r ■ 3.  For  r = Ij.  the  factor  (7.17)  exceeds  unity  in 
the  region  determined  by 


(7.19) 


.96 


and  is  less  than  unity  elsewhere*  In  view  of  the  remarks  made  about 
x,y,z)  and  since 


by  far  most  of  the  contribution  to  the  first  member  of  the  right  hand 
sides  of  (7.1U)  and  (7*15)  comes  from  the  integration  over  the  region 
determined  by  (7.19) • The  situation  is  similar  for  r «*  3*5  or  any  other 
value  of  r between  3 and  4.  Moreover,  in  the  neighborhood  of  the  mean 
the  contribution  to  the  first  member  of  the  right  hand  side  of  (7.1?)  is 
from  18  to  32  times  the  contribution  to  the  first  member  of  the  ri^ht 
hand  side  of  (7.1U)  for  values  of  r between  3 and  U<*  Consequently,  in 
all  practical  cases  Var  jj  S ^ar 

Since  the  computing  times  and  memory  requirements  are  about  the 
same  order  of  magnitude  for  the  three  random  sampling  integrating  pro- 
cedures the  smallest  variance  for  the  same  number  of  sampling  points  is 
the  deciding  criterion,  especially  in  view  of  the  fact  that,  in  order  to 
achieve  the  same  accuracy  with  the  methods  with  higher  variances  for 
the  same  number  of  sampling  points,  more  points,  and  hence  more  computing 
time,  are  needed.  Therefore,  RSIP-II  is  recommended  and  was  actually 
run  on  QRDVAC. 

The  two  systematic  numerical  integrating  procedures  will  be  con- 
sidered next. 

In  SNIP-I,  using  volumes  6a* s,  1atB,  and  8a's  on  an  edge  and 

choosing  5^,  and  6^  points  at  the  centers  of  equal  volume  elements 
whose  sides  are  parallel  to  those  of  the  large  volume,  the  kill  proba- 
bilities were  found  to  be  in  fair  agreement;  the  variations  were 
approximately  the  same  as  the  root  mean  square  relative  errors  in  the 
random  sampling  procedures.  However,  the  agreement  among  the  blast 
kills  was  poor;  this  was  apparently  due  to  the  fact  that  the  points 
chosen  lie  in  planes  parallel  to  the  plane  of  the  wings  of  the  air- 
plane and  are  separated  by  quite  a few  feet  thus  giving  a poor  sample 
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of  points  in  the  blast  regions*  Tried  out  in  a few  cases  with 
analytically  integrable  integrands  as  well,  SNIP- I gave  poorer  results 
than  SNIP-n.  Therefore,  SNIP-I  is  not  recommended  over  either  the 
random  sampling  procedures  or  SNIP-II,  even  though  its  computing  time 
is  about  two-thirds  that  of  the  random  sampling  procedures  for  the  same 
number  of  points*  SNIP-I  is  also  slightly  longer  than  SNIP-II  in 
computing  time* 

In  SNIP-II  1^,  and  (?  points  were  chosen  spaced  as  indicated 
in  the  previous  section;  these  were  run  for  all  five  sets  of  a *» 

For  the  trials  using  points  a majority  of  the  results  were  within  the 
observed  standard  deviations  of  the  means  in  RSIP-II  and  no  results  were 

•a 

in  bad  agreement  with  those  of  RSIP-II,  The  results  of  the  trials  using 

and  6^  points  were  very  well  within  the  observed  standard  deviations  of 
the  means  in  RSIP-II  in  almost  all  cases.  Furthermore,  the  results  using 

and  6^  points  agreed  very  closely  with  each  other,  in  several  ex- 
tremely good  cases  agreeing  to  within  three  or  four  units  in  the  third 
significant  figure.  In  general,  the  differences  between  the  results 

using  5^  and  points  were  about  half,  or  less,  the  observed  standard 
deviations  of  the  means  in  RSIP-II  using  averages  of  200  trials.  These 
results  are  taken  to  indicate  that  the  probabilities  obtained  by  SNIP-II 
using  125  burst  points  are  better  than  those  obtained  by  RSIP-II  using 
200  points.  Furthermore,  the  computing  time  for  SNIP-II  is  about  1.1* 
minutes  for  61*  points  for  each  set  of  a's,  about  2.2  minutes  for  120 
points  for  each  set  of  a's,  and  about  5 minutes  for  2l6  points  for  each 
set  of  a's,  while  in  RSIP-II  about  1*.0  minutes  are  required  for  100 
points  for  each  set  of  a's.  The  values  quoted  were  actually  timed;  for 
the  smaller  a's  the  blast  effects  were  higher  and  less  time  was  needed 
because  the  loops  involving  the  determination  of  kills  by  destruction 
of  components  were  omitted  for  many  points  but  the  times  were  higher 
whenever  the  blast  effects  were  low  for  the  converse  reason.  The 
results  of  these  investigations  indicates  that  for  one-missile  en- 
counters SNIP-II  is  recommended  over  SNIP-I  and  the  random  sampling 
procedures. 

II. 8 Basic  Logical  Flow  Chart. 

The  procedures  discussed  in  the  previous  sections  are  parts  of 
the  basic  program  whose  logical  flow  chart  Is  approximately  that  given 
in  Fig.  3 and  which  we  now  describe  briefly. 

Let  n be  the  Index  of  a burst  point  and  M be  the  total  number 

(usually  IDO  for  the  Lotto  method)  to  be  used  in  the  computation  of  I. 

Let  j be  the  index  of  the  type  of  vital  component  and  J be  the  total 

number  of  types.  Let  i be  the  index  of  a component  of  which  there  will 

be  I.  of  type  j.  Thus,  C.  . is  the  i-th  component  of  the  j-th  type, 
j 1,  j 

i^£l y j 
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Initially  n is  set  equal  to  1*  Then  a burst  point  Pn  is  produced 
as  a triplet  of  numbers  (xn,  y , z^)  generated  by  either  "random11 

sampling  or  by  systematic  methods  as  described  ini  H.6,  11*7*  and  the 
Appendix,  Then  it  is  determined  whether  or  not  the  point  is  inside  at 
least  one  of  the  ellipsoids,  the  sum  of  which  determines  the  blast 
region.  If  the  point  is  within  the  blast  region,  then  in  the  register 
reserved  for  the  integrands,  f(x,y,z),  a 1 is  stored  for  the  probability 
of  kill  from  a burst  at  P j otherwise,  the  probability  of  kill  is  deter- 
mined from  the  probability  of  kill  due  to  fragnents  as  follows. 


We  begin  with  determining  the  probability  of  destroying  the  first 
component  of  the  first  type  by  fragments  from  a burst  at  Pnj  thus, 

initially,  the  index  of  the  type,  j,  is  set  equal  to  1 and  i is  also 
set  eqial  to  1,  Then,  using  the  formulas  of  II.]*,  it  is  determined 

whether  C.  . is  in  the  conical  spray  regions  (See  Fig.  1,).  If  it  is, 
-*-*3 

then  it  is  determined  whether  or  not  C,  . is  shielded  by  any  of  the  K 

J 

ellipsoids,  E.  . If  C.  . is  not  shielded  by  any  of  the  K ellipsoids 
“ J 

then,  using  vulnerability  data  and  the  distance  between  (1  and  Pn, 

one  computes  and  stores  the  probability  of  destroying  C.  . by  frag- 

J 

ments  from  P . New,  if  i<  I,,  i,e«,  not  all  components  of  type  j 

have  been  examined,  i is  replaced  by  i + 1 and  the  loop  (or  an  alter- 
nate one  to  be  described:  below)  is  repeated  for  C,  _ If,  on  the 

i J 

other  hand,  i *»  I.  then  one  determines  whether  or  not  the  number  of 

unshielded  vital  components  of  type  j is  sufficient  for  a fragment 
kill  from  P^,  If  it  is,  then  one  computes  from  combinatorial  formulas 

and  stores  the  probability  of  a kill  due  to  destruction  of  type  j 
components  by  fragments  fran  P^.  If  it  is  not,  then  this  probability 

is  zero  and  0 is  stored. 


If,  on  the  other  hand,  C.  .is  not  in  the  spray  region  or  is 

3 

shielded  by  one  of  the  K ellipsoids  then  a 0 is  stored  for  the  proba- 


if  i' 


zi> 


bility  of  destroying  with  fragnents  from  P^.  Now, 

the  component  C. . replaces  C.  . for  a tour  through  either  this  loop 
a J.,  j if  j 

or  the  above  described  loop  and,  if  i = 1^  and  j < J the  procedure  is 

repeated  with  j + 1 replacing  j until  j ■ J.  Then,  using  the  stored 
probabilities  of  obtaining  kills  by  destroying  sufficient  numbers  of 
each  of  the  different  types  of  components  in  standard  combinatorial 
formulas,  the  probability  of  a kill  due  to  fragments  from  Pn  is  computed 

and  stored  in  the  register  reserved  for  f(x.j,y,z). 


Thus,  at  this  stage  we  have  in  the  register  far  f(x,  y,  z) 
the  probability  of  a kill  from  a burst  at  P^,  which  is  either  1 if  it 

is  a blast  kill  or  less  if  it  is  a fragment  kill.  This  procedure  is 
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repeated  until  n = M,  each  time  accumulating  f(xn,  yn»  zn)  in  the  same 

register.  Then,  multiplying  this  sum  by  the  proper  constants  (Since  I 
is  a probability  integral  and  since  different  methods  of  integration  may 
be  used  some  normalization  is  necessary,),  we  obtain  I,  which  is  then 
printed  and  the  program  halts. 

Many  things  have  not  been  included  in  this  flow  diagram  for  the 
sake  of  brevity  and  because  the  special  demands  of  each  program  intro- 
duce details  that  vary  from  problem  to  problem,  Segnents  of  this  flow 
chart  have  been  amplified  in  Fig*  2.  Things  such  as  determining  and 
printing  separately  the  probabilities  of  kill  of  various  components  have 
been  omitted  since  it  is  not  difficult  to  see  how  they  could  be  in- 
corporated in  this  program  (and  actually  have  been  in  several  problems 
solved* ) 


III.  CONCLUDING  REMARKS 

HI, 9 General  Remarks  on  the  Methods  and  Generalization. 

It  takes  about  ten  days  using  the  Lotto  method  by  hand  to  produce 
the  firing  record  tables  and  final  probabilities  for  five  different  sets 
of  variances  in  the  trivariate  normal  distribution  for  the  burst  point. 
Machine  time  on  QRDVAC,  which  included  printing  time  and  time  of  input 
from  IBM  cards,  to  produce  with  greater  accuracy  than  in  the  Lotto 
method  the  same  data  exclusive  of  the  firing  record  tables  varied  from 
3.2  to  3.8  minutes  using  SNIP-II  with  125  points  and  from  about  5 to 
5.5  minutes  using  RSIP-II  with  ICO  points.  The  variations  within  each 
method  were  due  to  the  variations  in  the  number  of  times  the  program  had 
to  investigate  fragment  kills  (i.e.,  when  the  burst  point  was  not  in  the 
blast  region),  the  fragnent  kill  loop  in  the  program  (See  Fig.  3.)  being 
much  more  complicated  than  the  blast  kill  loop.  RSIP-I  (the  Lotto 
method  on  the  machine)  would  take  slightly  longer.  However,  this  does 
not  include  coding  and  code-checking  time,  which  took  a few  weeks. 
Nevertheless,  since  it  is  foreseen  that  the  problem  is  of  a recurring 
type,  this  time  is  considered  as  initial  overhead  as  the  original 
formulation  of  the  Lotto  method  and  construction  of  models  must  also 
have  been.  For  subsequent  re-runs  of  the  problem  with  different  initial 
data  it  is  not  necessary  to  repeat  the  coding;  only  a small  amount  of 
time  is  required. 


An  additional  advantage  in  accuracy  in  the  ORDVAC  solution  over 
the  Lotto  method  appears  in  the  consistency  of  the  decision  as  to  when 
the  burst  is  in  the  blast  region.  In  the  Lotto  method  the  operators 
determine  by  eye  whether  or  not  a burst  point  is  in  the  blast  region, 
the  doubtful  cases  to  be  resolved  by  the  use  of  mathematical  formulas; 
but  for  psychological  reasons  these  seldom  occur.  (Mr.  King  has 
suggested  that  the  placement  on  the  model  of  a wire  mesh  form  outlining 
the  blast  region  would  eliminate  this  indeterminacy  in  the  hand  method.) 


On  the  other  hand,  in  the  formulation  of  the  problem  for  CRDVAC  the 
boundary  surface  of  the  blast  region  is  represented  by  a mathematical 
surface  of  second  degree  such  that  simple  discriminations  on  the  sign  of 
a quadratic  form  definitely  determine  whether  or  not  the  burst  is  inside 
the  blast  region. 


The  firing  record  tables  in  the  Lotto  method  have  an  advantage  in 
that  a change  in  kill  criteria  does  not  mean  that  new  tables  must  be 
made.  With  the  already  made  tables  within  several  hours  new  probabilities 
can  be  computed.  Furthermore,  the  acquisition  of  large  numbers  of 
firing  record  tables  represent  large  amounts  of  data  which  can  be,  and 
have  actually  been,  used  to  get  the  probabilities  of  kills  in  engage- 
ments involving  several  missiles.  (Of  course),  this  assumed  that  the 
cumulative  damage  due  to  bursts  from  different  missiles  on  a component 
is  not  enough  to  destroy  the  component  unless  the  damage  inflicted  by 
at  least  one  of  the  missiles  on  the  component  is  sufficient  to  destroy 
it.  This  would  also  bei  an  assumption  of  most  practical  formulations.) 

The  former  advantage  can  be  offset  somewhat  in  the  CRDVAC  because  of  the 
short  running  times  and  of  the  few  hours  necessary  to  produce  a tape  or 
cards  modifying  the  input  kill  criteria.  On  the  other  hand,  no  formu- 
lation for  several  missiles  has  been  attempted  for  CRDVAC  and,  there- 
fore, no  practical  comparison  between  the  methods  of  RSEP-I  and 
RSIP-II  are  now  possible j but  it  seems  likely  that  for  a large  number 
of  missiles  RSIP-I  would  be  preferred  over  RSIP-II. 

Further  generalizations  of  seme  sort  in  the  problem  seam  to  be 
possible  in  the  formulation  for  the  CRDVAC  solution.  A somewhat  more 
general  angular  distribution  of  fragments  is  possible  as  well  as  are  a 
distribution  of  directions  of  approach  of  the  missile  and  some  con- 
sideration of  the  directional  effect  of  blast  from  an  exploding  moving 
missile.  Clearly,  the  constants  in  the  formulas  for  the  probabilities 
of  destruction  of  components  by  fragnents  are  also  possible  to  be  changed 
with  little  delay  in  coding  or  computing. 

It  should  be  remarked  that  recently  there  has  been  added  the 
logical  "and”  to  the  list  of  automatic  CRDVAC  orders.  This  permits 
the  rapid  use  of  the  mare  general  method  described  first  in  II.5,  based 
on  the  partial  ordering  (<  ),  for  determining  whether  a certain  com- 
bination of  components  destroyed  constitutes  a kill  or  not.  The  method 
as  originally  coded  not  using  this  order  required  a very  large  amoung 
of  time  simply  for  shifting  in  order  to  make  the  digit-by-digit  com- 
parisons necessary  in  determining  whether  or  not  the  relation  (tf ) is 
satisfied.  The  elimination  of  this  shifting  makes  the  time  required 
about  the  same  as  that  required  for  the  less  general  method  described 
in  II.5- 

The  authors  wish  to  acknowledge  gratitude  to  Mr.  F.  G.  King  for 
many  preliminary  discussions  acquainting  them  with  the  problem,  to 
Dr.  Saul  Goto  for  helpful  suggestions,  in  particular,  iri  the  geometry 
of  the  problem,  and  to  Mr.  Frank  Lerch  for  machine  information  which 
helped  the  authors  crystallize  the  mathematical  formulation  of  the 
problem. 
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APPEND H 


10,  The  Generation  of  the  Random  Numbers . 

The  so-called  “random*1  numbers  used  in  the  randan  sampling  pro- 
cedures are  not  random  in  the  strict  mathematical  sense  of  the  word,  but 
rather  they  possess  to  some  degree  only  some  of  the  properties  of  truly 
random  numbers.  Density,  frequency  of  occurence  of  certain  digit  com- 
binations in  certain  positions,  and  contingencies  were  seme  of  the 
properties  investigated  and  compared  with  the  theoretical  behavior  of 

o 

truly  random  numbers  by  means  of %~  tests.  More  properly  these  random 
numbers  which  are  generated  in  a completely  deterministic  manner  are 
sometimes  called  pseudo-random  numbers,  but  for  economy  of  expression 
in  this  report  they  will  still  be  called  random  numbers. 

The  following  method  of  generating  the  random  numbers  directly  in 
our  problems  is  similar  to  the  procedure  described  by  D.  H.  Lehmer 
(p.  lUi  f 3]  ) and  to  the  procedure  devised  by  Olga  Taussky-Todd  for 
SEAC.  It  produces  sequences  of  numbers  approximating  uniformly  dis- 
tributed random  numbers  in  (0,l)  and  has  the  advantages  of  requiring 
very  few  orders  and  of  producing  sequences  having  a very  long  period 

(2"?f)  and  satisfying  very  well  certain  so-called  random  number  tests, 

39 

Let  p be  an  arbitrary  odd  number  satisfying  1 5 - 1. 

Define 


Pn+1  s 513  pn  (mod  239),  n e 0,  1,  2,  ..., 

39 

and  such  that  1 1 p^S  2 - 1 for  every  n*  Then  define 

rn  = 2-39  pny  n - 0,  1,  2,  ...  . 

The  sequence  ^ rn^  is  the  desired  random  number  sequence  whose  distri- 
bution is  very  close  to  the  uniform  distribution  in  (0,l).  On  CRDVAC 
this  is  simply  achieved  by  multiplying  rR  by  $ J using  double  precision, 
i.e.,  using  two  registers  (78  binary  places)  for  the  product,  and  in 
such  a manner  that  the  integral  part  of  $ J rn  falls  in  one  register 
and  the  fractional  part,  which  is  rn+^,  falls  in  the  other  and  is 
ready  to  be  used  in  the  problem  as  well  as  to  generate  rn+2* 

39 

The  modulus  2J  was  chosen  because  the  CEDVAC  register  has  39 
binary  digits  and  the  obtaining  of  the  remainder  of  the  division  of 
13  39 

$ by  2J  is  achieved  immediately  by  omitting  the  first  register  in 

the  double  precision  division.  The  multiplier  £^3  was  chosen  because 

39  13 

it  does  not  exceed  2J7  and  thus  $ 7 pn  fits  entirely  in  the  two 
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13 

registers  reserved  for  it;  on  the  other  hand,  $ has  a reasonable 
selection  of  digits  and  is  large  enough  for  and  pn  not  to  appear 

correlated.  Furthermore,  in  order  to  insure  the  longest  possible 

period  with  this  sort  of  scheme  the  multiplier  should  be  congruent  to 

13 

1 with  respect  to  the  modulus  1*,  which  requirement  is  satisfied  by  5 . 

It  is  perhaps  worthwhile  to  remark  that  this  process  is  easily 
generalized  for  any  high-speed  digital  computing  machine.  For  a machine 

whose  registers  contain  n digits  to  the  base  p,  the  modulus  should  be  pn 

instead  of  2^;  the  multiplier  used  in  place  of  should  be  large  tut 

less  than  pn  and  if  p = 2 the  multiplier  should  be  congruent  to  1 with 
respect  to  the  modulus  1*  in  order  to  insure  the  longest  possible  period 

(P  “ ),  If  p ^ 2,  a special  analysis  for  each  machine  is  necessary  to 
determine  to  which  residue  classes  the  multiplier  should  belong* 

Using  elementary  number- theoretic  methods,  one  can  show  that  this 

procedure  on  ORDVAC  will  produce  a succession  of  exactly  2^  different 
odd  numbers,  pn  and  will  then  continue  to  repeat  the  se  (pence  over 

again.  This  implies  that  exactly  half  of  all  the  odd  numbers  in  the 

interval  1^  pn^  2 - 1 will  appear  in  each  period;  the  particular 

set  of  odd  numbers  appearing  will  depend  on  the  particular  choice  of  pQ, 
Furthermore,  one  can  easily  show  that,  if,  for  same  n,  pn  = p,  then 

there  exists  no  k such  that  pn+j^  = p + 2 in  the  same  sequence:  thus, 

the  two  sets  of  odd  numbers  which  may  be  produced  interlace. 

Using  the  iterative  procedure  pn+^  e Pn  (mod  2^),  pQ  - 1,  to 

produce  pseudo-random  numbers  on  the  SEAC,  the  National  Bureau  of 
Standards  made  seme  fairly  extensive  and  exhaustive  tests  whose  re- 
sults indicated  very  good  agreement  with  what  cculd  be  expected  from 
truly  random  numbers.  Since  our  method  is  very  similar  and  their  re- 
sults were  so  good,  not  so  extensive  tests  were  performed  on  the 
sequences  produced  on  ORDVAC. 

With  rQ  = 1 - and  rQ  ■ 0.51*78126193  two  seepences  of  U096 

numbers  were  produced.  The  number  of  zeros  in  each  of  the  following 
places  of  the  binary  representations  of  the  numbers  of  the  sequences : 
2nd,  3rd,  5 th,  7th,  11  th,  18th,  and  21*th,  was  counted  for  each  sequence. 
2 2 

Values  of  V ■ 2,87  and  ■»  7,1*8  respectively  were  obtained.  These 

O Op 

values  are  exceeded  by  a a variable  with  seven  degrees  of  freedom  with 
probabilities  of  0.82  and  0,39  respectively.  The  number  of  occurrences 
of  10  in  the  l*th  and  5th  and  in  the  l£th  and  16th  places  of  the  binary 
representations  of  the  numbers  in  each  sequence  was  counted,  giving 
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19 

1067  and  1057  for  the  sequence  beginning  with  r = 1 - 2J7  and  1060 

and  103?  for  the  sequence  beginning  with  rQ  *»  0 .5478126193 « The 

probability  of  random  fluctuations  from  1024>  the  expected  number  of 
such  occurrences,  exceeding  43/  the  largest  of  the  above  deviations,  is 
0.f>0  using  the  central  limit  theorem  for  4096  trials.  The  number  of 
occurrences  of  1011  in  the  7th  through  10th  places  of  the  binary  repre- 
sentation of  each  number  was  counted  for  each  sequence,  giving  deviations 
of  4 and  18,  respectively;  the  probabilities  that  such  deviations  would 
be  exceeded  in  truly  random  sequences  are  0.9?  and  0.78  respectively. 

An  investigation  of  the  density  of  the  sequence  of  4096  numbers 

produced  from  r ■ 1 - 2“-^  was  made.  Dividing  in  the  interval  into 

200  equal  sub-intervals  a value  of  of  199*4L  was  obtained;  this  is 

exceeded  by  a ^ variable  with  199  degrees  of  freedom  with  a probability 

of  0.46.  A different  division  of  the  same  sequence  into  10  equal  sub- 

2 2 
intervals  gaveX0  ■ 6,95,  which  is  exceeded  by  aX  variable  with  9 

degrees  of  freedom  with  a probability  of  0.58*  With  the  same  division 
into  10  equal  sub-intervals  a contingency  table  for  rn  and  rn+1  and  for 

rn  and  rn+^  producing^  = 88.03  and^  “ 92.33*  respectively,  which  are 

exceeded  by  X?  variable  with  99  degrees  of  freedom  with  probabilities 
of  0*78  and  O.67  respectively.  On  the  basis  of  these  tests  and  the 
modest  accuracy  requirements  in  our.  problem  this  procedure  was  accepted 
as  a method  of  generating  sequences  of ; pseudo-random  numbers  approxi- 
mating sequences  of  uniformly  distributed  random  numbers  in  (0,l). 

Furthermore,  a X < test  was  made  by  F.  G.  King  on  the  values  obtained  far 
the  probability  of  a blast  kill  in  his  requested  $0  runs  of  the  problem 
using  50  independently  randomly  chosen  ro*s  and  he  found  the  numbers 

acceptable  for  the  purposes  of  the.  problem. 

To  produce  the  pseudo-random  normally  distributed  numbers  required 
in  the  random  sampling  procedures  used  in  the  problem  an  elementary  and 
well-known  device  was  used.  If  Y is  a uniformly  distributed  random 
variable  and  if  F(x)  is  a strictly  monotone  increasing  cumulative  dis- 
tribution function  whose  inverse  is  F^Cx),  then  X » F”^(Y)  is  a random 
variable  with  cumulative  distribution  function.  F(x).  Thus,  for  our 
purposes  a simple  approximation  to  the  inverse  of  the  normal  distri- 
bution was  used;  it  is  based  on  the  f omnia  approximately  inverting 


given  on  sheet  67  of  Form  (15 )s  £l]. 

If  Y is  a pseudo-random  number  approximately  a uniformly-distri- 
buted number  in  (0,l)  and 

Z - V-  2 log^  \ (1  - |l  - 2 Y|) 
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then 


X = m + 


Z 


a + a_  Z 

2 J: a sign  (Y  - .$) 

1 + Z + b2  Zc  * 


is  a pseudo-random  normal  variable  with  mean  m and  variance  a . The 
constants  aw  given  by 


aQ  - 2.307*3  \ = 0,99229 

ax  - 0,27061  b2  « O.OUiSl 


This  procedure  of  converting  uniformly  distributed  random  numbers  to 
normally  distributed  random  numbers  requires  about  80  words  of  Btorage 
and  takes  about  one- third  of  the  computation  time*  F.  G.  King  has 
pointed  out  that  a sum  of  eight,  or  possibly  even  four  with  a correction 
factor,  uniformly  distributed  random  numbers  as  would  be  a sufficiently 
close  approximation  to  a normally  distributed  random  number  for  the 
purposes  of  his  problem.  It,  furthermore,  would  require  less  storage 
and,  although  more  random  numbers  need  to  be  generated,  would  require 
less  computing  time.  Thus,  under  strait  circumstances'  it  would  be  ad- 
visable to  use  this  scheme  rather  than  the  one  now  coded.  One  minor 
consideration  is  necessaryj  more  extensive  testing  of  the  pseudo-randcm 
sequences  is  needed,  there  being  several  times  more  random  numbers  re- 
quired. 


This  latter  procedure  for  obtaining  normally  distributed  numbers  from 
uniformly  distributed  ones  was  coded  in  a later  revision  of  the  problem. 

It  used  a sum  of  5 or  6 uniformly  distributed  numbers  properly  normalized 


and  with  a small  correction  to  improve  the  ta, 
time  saving  in  the  running  of  the  problem 
amounted  to  about  30 # to  to#  of  the  previous 
time.  This  is  due  to  the  elimination 
of  the  square  root  and  logarithm  routines 
needed  In  the  former  procedure. 


the  distribution.  The 
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